{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Data poisoning attack\n",
    "\n",
    "In this notebook, we use a convex optimization layer to perform a *data poisoning attack*; i.e., we show how to perturb the data used to train a logistic regression classifier so as to maximally increase the test loss. This example is also presented in section 6.1 of the paper *Differentiable convex optimization layers*."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "import cvxpy as cp\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import torch\n",
    "\n",
    "from cvxpylayers.torch import CvxpyLayer"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We are given training data $(x_i, y_i)_{i=1}^{N}$,\n",
    "where $x_i\\in\\mathbf{R}^n$ are feature vectors and $y_i\\in\\{0,1\\}$ are the labels.\n",
    "Suppose we fit a model for this classification problem by solving\n",
    "\\begin{equation}\n",
    "\\begin{array}{ll}\n",
    "\\mbox{minimize} & \\frac{1}{N}\\sum_{i=1}^N \\ell(\\theta; x_i, y_i) + r(\\theta),\n",
    "\\end{array}\n",
    "\\label{eq:trainlinear}\n",
    "\\end{equation}\n",
    "where the loss function $\\ell(\\theta; x_i, y_i)$ is convex in $\\theta \\in \\mathbf{R}^n$ and $r(\\theta)$ is a convex\n",
    "regularizer. We hope that the test loss $\\mathcal{L}^{\\mathrm{test}}(\\theta) =\n",
    "\\frac{1}{M}\\sum_{i=1}^M \\ell(\\theta; \\tilde x_i, \\tilde y_i)$ is small, where\n",
    "$(\\tilde x_i, \\tilde y_i)_{i=1}^{M}$ is our test set. In this  example, we use the logistic loss\n",
    "\n",
    "\\begin{equation}\n",
    "\\ell(\\theta; x_i, y_i) = \\log(1 + \\exp(\\beta^Tx_i + b)) - y_i(\\beta^Tx_i + b)\n",
    "\\end{equation}\n",
    "\n",
    "with elastic net regularization\n",
    "\n",
    "\\begin{equation}\n",
    "r(\\theta) = 0.1\\|\\beta\\|_1 + 0.1\\|\\beta\\|_2^2.\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.datasets import make_blobs\n",
    "from sklearn.model_selection import train_test_split\n",
    "\n",
    "torch.manual_seed(0)\n",
    "np.random.seed(0)\n",
    "n = 2\n",
    "N = 60\n",
    "X, y = make_blobs(N, n, centers=np.array([[2, 2], [-2, -2]]), cluster_std=3)\n",
    "Xtrain, Xtest, ytrain, ytest = train_test_split(X, y, test_size=.5)\n",
    "\n",
    "Xtrain, Xtest, ytrain, ytest = map(\n",
    "    torch.from_numpy, [Xtrain, Xtest, ytrain, ytest])\n",
    "Xtrain.requires_grad_(True)\n",
    "m = Xtrain.shape[0]\n",
    "\n",
    "a = cp.Variable((n, 1))\n",
    "b = cp.Variable((1, 1))\n",
    "X = cp.Parameter((m, n))\n",
    "Y = ytrain.numpy()[:, np.newaxis]\n",
    "\n",
    "log_likelihood = (1. / m) * cp.sum(\n",
    "    cp.multiply(Y, X @ a + b) - cp.logistic(X @ a + b)\n",
    ")\n",
    "regularization = - 0.1 * cp.norm(a, 1) - 0.1 * cp.sum_squares(a)\n",
    "prob = cp.Problem(cp.Maximize(log_likelihood + regularization))\n",
    "fit_logreg = CvxpyLayer(prob, [X], [a, b])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Assume that our training data is subject to a data poisoning attack,\n",
    "before it is supplied to us. The adversary has full knowledge of our modeling\n",
    "choice, meaning that they know the form of the optimization problem above, and seeks\n",
    "to perturb the data to maximally increase our loss on the test\n",
    "set, to which they also have access. The adversary is permitted to apply an\n",
    "additive perturbation $\\delta_i \\in \\mathbf{R}^n$ to each of the training points $x_i$,\n",
    "with the perturbations satisfying $\\|\\delta_i\\|_\\infty \\leq 0.01$.\n",
    "\n",
    "Let $\\theta^\\star$ be optimal.\n",
    "The gradient of\n",
    "the test loss with respect to a training data point, $\\nabla_{x_i}\n",
    "\\mathcal{L}^{\\mathrm{test}}(\\theta^\\star)$, gives the direction\n",
    "in which the point should be moved to achieve the greatest\n",
    "increase in test loss. Hence, one reasonable adversarial policy is to set $x_i\n",
    ":= x_i +\n",
    ".01\\mathrm{sign}(\\nabla_{x_i}\\mathcal{L}^{\\mathrm{test}}(\\theta^\\star))$. The\n",
    "quantity $0.01\\sum_{i=1}^N \\|\\nabla_{x_i}\n",
    "\\mathcal{L}^{\\mathrm{test}}(\\theta^\\star)\\|_1$ is the predicted increase in\n",
    "our test loss due to the poisoning."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.linear_model import LogisticRegression\n",
    "\n",
    "a_tch, b_tch = fit_logreg(Xtrain)\n",
    "loss = 300 * torch.nn.BCEWithLogitsLoss()((Xtest @ a_tch + b_tch).squeeze(), ytest*1.0)\n",
    "loss.backward()\n",
    "Xtrain_grad = Xtrain.grad"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Below, we plot the gradient of the test loss with respect to the training data points. The blue and orange points are training data, belonging to different classes. The red line is the hyperplane learned by fitting the the model, while the blue line is the hyperplane that minimizes the test loss. The gradients are visualized as black lines, attached to the data points. Moving the points in the gradient directions torques the learned hyperplane away from the optimal hyperplane for the test set."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAD8CAYAAABjAo9vAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3dd3hU1dbH8e9OSAihCyhNuhSlSQcVpFcRBQkiClhoCogFVOyvnauiIgIqipWuiCJSBPVKUUCkVxtVQC69J+v9YyWCGEiZcmYm6/M8echMZs5ZRPxlZ5991nYigjHGmMgR5XUBxhhj/MuC3RhjIowFuzHGRBgLdmOMiTAW7MYYE2Es2I0xJsL4Jdidc4Occ6udc6uccx875+L8cVxjjDEZ53OwO+eKAQOAWiJSGYgGuvh6XGOMMZnjr6mYbEAO51w2IB7Y7qfjGmOMyaBsvh5ARLY55/4D/AEcBWaJyKyzX+ec6wX0AsiZM2fNihUr+nrqiCICmzfD/v1QogQUKhTkAk4ehd3roFBFiMkR5JMbY9Jj6dKle0QkzXRwvrYUcM7lB6YACcA+YBIwWUQ+ONd7atWqJUuWLPHpvJHo+HHo2BG++AJGj4ZevYJw0l1rYcrtJB3cxamDu4jNexHkLAQd34ILKwWhAGNMejnnlopIrbRe54+pmGbAryKyW0ROAlOBBn44bpaTPTtMmQJt2kDv3vDmm0E4ab6SbKEY0UM28ebSE3B4DxSqAPlKBuHkxphA8Eew/wHUc87FO+cc0BRY64fjZkkp4d66tY7Y33orsOdbsW4TJfpNAqBXnZwgidBwMMTGB/bExpiA8TnYRWQxMBlYBqxMPuYYX4+blcXFwdSp0KoV3HEHjB0bmPPMmTOHatWqER0dRdIbVxHTfSoUqQ67VgfmhMaYoPB5jj0zbI49fY4dgw4dYNYsePtt6NnTf8d+//33ueWWW6hSpQorVqzw34GNCZCTJ0+ydetWjh075nUpARcXF0fx4sWJiYn5x/PpnWP3eVWMCZy4OPj0U7j2WrjtNnAOevTw/bhPP/00Dz/8MJ06dWLSpEm+H9CYINi6dSu5c+emVKlS6KxvZBIR/vrrL7Zu3Urp0qUzdQxrKRDiUsK9WTO49VYYN86349166608/PDDDB482ELdhJVjx45RoECBiA51AOccBQoU8Ok3Exuxh4EcOWDaNB259+ypI/dbbsn4cd5++23eeecdRo4cSd++ff1fqDEBFumhnsLXv6cFe5hICff27XU6xjm4+eaMHSMhIYF69epx2WWXBaRGY0xosKmYMJIS7k2aQPfu8ME5bwFLXa5cuSzUjfHBvn37GDlyZIbf16ZNG/bt2xeAilJnwR5m4uPhs8+gcWMN9w8/9LoiY7KOcwV7YmLied83Y8YM8uXLF6iy/sWmYsJQfDxMnw7t2ulcu3PQtavXVRkT+R544AE2b95M9erViYmJIVeuXBQpUoTly5ezZs0aOnTowJYtWzh27BgDBw6kV3JfkFKlSrFkyRIOHTpE69atufLKK1mwYAHFihVj2rRp5Mjh3/5MFuxh6sxwv/lmDfcbb/S6KmOC6Oqr//1c587Qrx8cOaK9Oc7Wo4d+7NkDnTr982vz56d5yueee45Vq1axfPly5s+fT9u2bVm1atXfyxLHjh3LBRdcwNGjR6lduzYdO3akQIEC/zjGxo0b+fjjj3nzzTfp3LkzU6ZMoVu3bun6K6eXTcWEsZw54fPP4aqroFs3mDDB64qMyVrq1Knzj7Xmr776KtWqVaNevXps2bKFjRs3/us9pUuXpnr16gDUrFmT3377ze912Yg9zOXMqd0g27SBm27SkXvnzl5XZUwQnG+EHR9//q8XLJiuEXpacubMeUY585kzZw4LFy4kPj6eq6++OtW16NmzZ//78+joaI4ePepzHWezEXsESAn3Bg10rn3iRK8rMiYy5c6dm4MHD6b6tf3795M/f37i4+NZt24dixYtCnJ1p9mIPULkygUzZmhXyK5dISrq31OIxhjfFChQgCuuuILKlSuTI0cOLrroor+/1qpVK0aNGkXVqlWpUKEC9erV86xOawIWYQ4e1HBftEjn3Dt29LoiY/xj7dq1VKqUdTZ/Se3vG8yNNjJu/Xr46CNIY+2nybjcueHLL6FuXejSRdv/GmOyFm+CPTFRr/RVqaITwklJnpQRqVLCvXZtSEiATz7xuiJjTDB5E+yXXqqB7pwmT7VqsHevJ6VEqjx5YOZMqFVLV8l8+qnXFRljgsW7VTE33AArVuiUTO3akD+/Pr96NXgw7x+JUsK9Zk39dk+b5nVFxphg8Ha5Y3S03i45dqyO3nfs0CFm7dq6fs8C3md588JXX50O988+87oiY0yg+SXYnXP5nHOTnXPrnHNrnXP1M3WgQoVg5EidlmnXDurV01SygPdJSrhffrkugZw+3euKjDGB5K8R+yvATBGpCFQD1mbqKNmy6U4S69fDmDGwc6feUvnrr34qM+tKCffq1XUJ5Oefe12RMeEns217AYYPH86RI0f8XFHqfA5251weoCHwNoCInBAR3xoPx8TAHXfAxo2aRmXK6PP/93/wzTc+Vpx15cunG2NXq6bh/sUXXldkTHjJMsEOlAF2A+84535yzr3lnMuZ1pvSJTZWN/sE2L8fRo3Sjm5Nm8J//+uXU2Q1KeFetSpcf73erWqMSZ8z2/bef//9DBs2jNq1a1O1alUee+wxAA4fPkzbtm2pVq0alStXZsKECbz66qts376dxo0b07hx44DX6Y+WAtmAGkB/EVnsnHsFeAB45MwXOed6Ab0ASpQokfGz5M0LmzbB6NHw3HPa0rB5c3jjDShb1ue/RFaSP7+Ge/PmcN11uhSydWuvqzIm/e6+G5Yv9+8xq1eH4cPP/5oz2/bOmjWLyZMn88MPPyAitG/fnm+//Zbdu3dTtGhRvkj+lXj//v3kzZuXl156iXnz5lGwYEH/Fp4Kf4zYtwJbRWRx8uPJaND/g4iMEZFaIlKrUKFCmTtTjhz6X/SXX2DYMA36PHn0a0H6FSdS5M8Ps2dD5coa7jNnel2RMeFl1qxZzJo1i8svv5waNWqwbt06Nm7cSJUqVZgzZw5Dhgzhu+++I2/evEGvzecRu4jsdM5tcc5VEJH1QFNgje+lnUd8PNx3HwwapEsmRXSvuIsugiee0OUfJk0p4d6sGXTooOvcW7b0uipj0pbWyDoYRIQHH3yQ3r17/+trS5cuZcaMGTz44IO0aNGCRx99NKi1+WtVTH/gQ+fcCqA68Iyfjnt+0dH656lT0L49fPcd1Kihk8crVgSlhHB3wQUwZw5UqgTXXqtTNMaY1J3Ztrdly5aMHTuWQ4cOAbBt2zZ27drF9u3biY+Pp1u3btx3330sW7bsX+8NNL8Eu4gsT55mqSoiHUTkf/44brrFxMDQofDbb/D44zB3ri79sGUf6ZIS7hUrarjPnu11RcaEpjPb9s6ePZuuXbtSv359qlSpQqdOnTh48CArV66kTp06VK9enaeffpqHH34YgF69etG6deugXDyNzLa9//ufrqAZNAji4uDrr6FoUU0uc0579uiCow0b9CamlAVJxoQCa9sb6m17Ay1/fnjwQQ11ERgwAC67THd9TmUPQqMKFtRfdi65BK65Rj83xoSfyAz2MzmnI/Z77oEpU3QyuWdPu5v1HM4O96+/9roiY0xGRX6wA1x4oS6P/OUXHb2PH28XV8+jUCEN97JltWWPhbsJFV5MHXvB179n1gj2FIULw0sv6Wi9fXt97qmnoG9f2LLF29pCTEq4lymj4e6HDd2N8UlcXBx//fVXxIe7iPDXX38RFxeX6WNkzc2sCxc+/fmhQ/D229o6uFcvnZsvWtS72kLIhRfqaL1xY2jbVtsPNGrkdVUmqypevDhbt25l9+7dXpcScHFxcRQvXjzT74/MVTEZ9fvv8MwzGu7R0fD663DbbV5XFTL+/BOaNNHVpBbuxngna6+KyaiSJbUHzYYNuhdr9er6/I4dsGuXt7WFgIsu0pF7yZLaRfnbb72uyBhzPhbsZypdWqdlatbUxw89pM8NGaKLvLOwlHAvUULD/bvvvK7IGHMuFuzn88ADeivmsGEa8A8/nKU33S5cGObNg4sv1m6Q1jnZmNBkwX4+FSroZturVukw9emntclYFla4sI7cixe3cDcmVFmwp8ell8KECbr2fcgQfW7hQt3R6cABb2vzQJEiOnIvWlTDfcECrysyxpzJgj0jqlQ5vRRy1ix49FGdonn2WV02mYWkhHuRItCqlf6cM8aEBgv2zHrsMfjxR6hf//RF1lGjvK4qqIoW1XAvXFj7uC9a5HVFxhiwYPdNrVrw+eeaaDVqaFdJgMREOHrU29qCpFgxDfeLLtJwX7w47fcYYwLLgt0f6taFr746Pf8+caI2WnntNTh2zNvagiAl3AsVghYt4IcfvK7ImKzNgt2fopK/nWXKQPny2nCsXDkYORKOH/e2tgArXlzDvWBBDfcff/S6ImOyLgv2QKhbV1Nu7lwoVQruvFObrUS4iy/WZmEFCkDz5hbuxnjFgj1QnNMGK999pyto7r9fnz98GMaNg5Mnva0vQC6+WH+mXXCBhnsotQQyJqvwW7A756Kdcz855z731zEjgnOacC1b6uPx46FHD93w4733dCPuCFOihIZ7/vz6V1+61OuKjMla/DliHwis9ePxItOtt8K0aZA7N3Tvrlv2ffSRbuEXQUqW1GmZfPl071QLd2OCxy/B7pwrDrQF3vLH8SKac7rJx9KlMHUqZM8OY8Z4XVVAlCypI/e8eXXkvmyZ1xUZkzX4a8Q+HBgMJJ3rBc65Xs65Jc65JVmhUX6aoqLguutg+XKYNEkDf/t2qFNHAz/pnN/KsFKqlI7cc+fWkftPP3ldkTGRz+dgd861A3aJyHl/2RaRMSJSS0RqFSpUyNfTRo6oKF0ADrBtm/ae6dhRb3iaNi0ipmjODvfly72uyJjI5o8R+xVAe+fcb8B4oIlz7gM/HDfrqV0bVq+G99/X1TMdOkC9enDihNeV+ax0aZ2WyZkTmja1cDcmkHwOdhF5UESKi0gpoAvwtYh087myrCo6Grp1g7Vrdau+Jk0gNla/9tNPYT2CL1NGwz0+XsP955+9rsiYyGTr2ENVtmzQs6d2jgRNwRo14Mor9canMA34smV1WiYl3Fes8LoiYyKPX4NdROaLSDt/HtMkq1QJ3ngD/vhDJ6qvvhq++cbrqjKlbFkducfFabivXOl1RcZEFhuxh4vYWOjTBzZt0uZiGzdqm4J9+7yuLF0SRi8kYfTppu3lymm4Z8+us02rVnlYnDERxoI93GTPDnfdBZs3a6uCfPl0WmbgwJBtiJ4weiGLf93L4l/3/iPgL7lEwz02VsN99WqPCw01IrDms7CddjPesWAPVzlyQIMG+vnvv8PHH+umH23bhlyDlv/9sZ7fn2/H/kWT//W1lHDPlk3Dfc0aP5wwUgJx5wqYeDPszORcVaR8H0yGWbBHglKl4Jdf4LnndNReu7be3bpzp6dliQgdO3Zk9jM9AWhy/c1M6F2fCb3r/+N15cvrBdXoaGjc2A/h7msgem3XWnjjCn4e1h5w8GEneOMKfT4jwv37YDLNgj1S5MqlG338+is89ZSO4vPl068dPBj0cjZt2kRUVBRTp05lxIgRdB61gGzZc5zz9eXL68g9KkpH7msz03UoORDXvNSev45I5gPRa/lKMmD6/6j+8h8cO5UEh/dAoQqQr2T63p/8feDDG/DpB0MKG/mHHQv2SJMnDwwdqncAxcXpBh9VqkDnzkGbxB4wYACXXHIJAHv27OHOO+9MdaR+tgoVNNyd05H7unUZPHG+koxalsRl//mD1bsTMx6IIWLEmLG89uVaXmyRnbjs2UESoeFgiI1P3wHyleRkvnK4+zcy/7eTvn8fbOQfdizYI5Vz+uepU3DLLTBzpgZ8166ZSMz02blzJ845XnvtNYYOHYqIUKBAgQwdo2JF+Ppr/Tyj4f7C8BH0fWsRt14eQ8My8RkPxBAwb948+vfvzx2dmnPP9XWh2xQoUh12neeH8tkj6th4Yru+B0C9Ejky/31IHvmffK8TR04Svr8BZUUiEvSPmjVrigmyPXtEHnhAJGdOkagokSVL/Hr4Z599VgAB5I8//vD5eKtXi1x4oUjhwiLr1qX9+iFDhggg9/W4VmRUQ5FfvtE/V072uZZgmjt3rvTv3z9jb9q+XOSxPCLbfxYRkZMnTwogKx+p6dv34fhhkYk9BJBHGsaKPJ5fZFJPfd54Algi6chYC/asZtcukRdfFElK0seffCKyeXOmD7dv376/A71nz55+KlKlhHuRIucP9969ewsgTz31lF/PH/L+XCMysoE83aqAdLo0RmTYJSIjG+jzfpK4Y5WULxAlW+/Lrz88/lzrt2ObjEtvsNtUTFZTqBDcc49O1Zw4Ab176+T2HXfoBdcMeOedd8iXfIF29erVjB071q+lXnqpTsucOqXTMhs2/Ps1Xbp0YfTo0X9P/2Qp+Uryza48DJ35FyXzuoBcU4jas471T19JsX6fpj0lZEJHetLf3x82Yg8h27aJ9O8vEhsrEhMj0qePyNat533L0aNHJS4uTgBp2bKlJKWM/gNk5UqRQoVEihYV2bDh9PPNmzcXQN59992Anj9UnTp1SgApmdeJPFnQRtRZADZiN+lStCi8+qreyXr77dpR8s8/z/ny6dOnkyNHDo4dO8aCBQuYOXMmLuVCbYBUrqx9z06c0JH7pk1Qs2ZNZs+ezdSpU+nevXtAzx+qoqOjmfvWY/z63JXpu8hqsgwnHqxNrVWrliwJsbsjTbI9e6BgQf28b19tYfDAAyQWKkSlSpXYuHEjlSpVYuXKlURHRwe1tBUrdI37vn07SEy8iq++GkmLFi2CWkOGiMDa6VDpmtOrlIzxgXNuqYjUSut1NmI3/5QS6iKQmAgjRpBYqhQvZ8vGvo0b+eKLL1izZk3QQx2gShXh4MHaJCbGcOGFqylbNoRDHWz9t/GMBbtJnXMwZgz3tm3LB8ePMwj4Mz6eNh7u5lSjRg1OnFjC+PF7SEzMTuPG2kkh5CSv/7653VVcO/6Irf82QWfBbs4rV/Xq8O67RK9bh+vYEWol/xa4eTPs3RvUWt544w02btxIQkJF5s7V3QOvvjoEwz1fSe77cj8fLDvItRWyhe0dsCZ82Ry7yZxmzeDHH+Huu2HQoNN9aYJo+XKdc8+dW5uIlS4d9BJSNWzYMAYPHszzzbIzuGFuSDwB/RbDhRW9Ls2EOZtjN4H18ssa7k8+qYn61FNw4EBQS6heXVfLHDyoq2V++y2op0/Vu+++y+DBgxl0yzUM7pTOlgDG+Ft61kSe7wO4GJgHrAVWAwPTeo+tY48gP/0kcu21ehPzCy94UsLSpSL584uULCny22+elCAiItOmTRNAbrzxRu+KMBGNdK5j93kqxjlXBCgiIsucc7mBpUAHETlnV22biolAS5fqHay5csGnn+rWff36Qc6cQTn9smW6f2q+fDotU9KD6WznHI0aNWL+/PnBP7nJEoI2FSMiO0RkWfLnB9GRezFfj2vCTM2aGuoAs2fD4MFQpoxO2Rw9GvDT16gBc+boFrCNG+ue38F29OhRC3UTEvw6x+6cKwVcDixO5Wu9nHNLnHNLdu/e7c/TmlDz+uvw/fdQtar2pSlbVrfuC7CaNfVnyt69ulpmy5aAn/If4uLigntCY87Bb8HunMsFTAHuFpF/XUUTkTEiUktEahUqVMhfpzWhqkEDTdlvvtHtkZKS9Pljx3TzjwCpVUtP+9df3oS7MaHAL8HunItBQ/1DEZnqj2OaCNGwoU56d+2qj195RXewHjNGm78EQO3aGu579ui0zNatATmNMSHL52B32gHqbWCtiLzke0kmIqX0SqlTRxuPpbQLHjsWTp70++nq1IFZs2DXLh25b9vm91MYE7L8MWK/ArgZaOKcW5780cYPxzWRqHFjWLgQZszQvjS33aYfAVC3roW7yZr8sSrmvyLiRKSqiFRP/pjhj+JMhHIOWreGH36Azz6D/v31+W3b4MMPtfmYn9SrB199pZ2IGzeOwHCXs/Y7NYETRt9ru/PUnF8g/zE7B9dco5PiAO+8A9266abbEyacvuDqo/r1dS/vHTu0BcH27X45bGiwDpLBE0bfawt2c37B/Mf80EMwaZIGfpcuUK2a3uzkBw0aaLhv364j9x07/HJY7yR3kGzduAGjlpywDpKBlPy9/mvM9eR4+gAf3tsy5L/XFuwmdbvWcuClurii1Vm9Kyk4wREVBZ066Y4aH3+sF1WnnrHIysffGq64Ar78Uqdjwj7c85Xk+z15mbn+CJUKRlkHyUDKVxIKlueXbbtxQMXch0P+e23BblJ1JLYQee/9AYCyFwRmo+Rzio7WEfvq1bptH8BPP+lSly++8Cngr7xSR+5bt+q0zM6dfqo52GLj6TxyKblioVHZeJBEaDgYYuO9rizyxMZDo8HULuo48mhBahZxIf+9tmA3/3L8+HFy5tebyA4+mJu47Nm9CY7o6NPtgPfu1buO2rU7fUU0kwF/5ZU6ct+yRUfu4Rruqya/wO4Xbb/ToNi1Rr/HYfK9tn7s5h9OnjxJbGwsAPsWjCPvirehxf/BrEfgigFQuaOXxcG4cdoi+PffoWVLTehM7if67bfQpg2UKAHz5sFFF/m5XmP8LL1NwCzYzd8SExPJli0bALt376Zgyv6noebECb2x6dgx3ehDRJdO1q2b4UN9842Ge6lS8PXXFu4mtNlGG36WMHohCaMXel1GwIjI36G+ffv20A11gNhY6NNHQx30Zqd69bRv73//m6FDNWqk0/a//aZz7rt2+b9cY4LNgj0dWj/xMVMf6syizbsjMuBFhKgo/afw+++/U6RIEY8ryqAmTWD4cL3YetVV0KKF3t2aTldfreH+6696KGs+mrb33nuPwUOGeF2GOQcL9jQ8+eSTzHy8K6f2bs30XG6oSwn1DRs2UKJECY+ryYQcOWDgQN3V+j//0c1QExIy1IPm6qvh88/1EGmGexjdgRgoK1asYNgLL0TcICdSZPO6gFB14sQJsmfPDsDAgQPZUSkBgAm963tZlt/lzp0bgFWrVnHJJZd4XI2P4uPh3nu1wdjGjRAToy2C+/WDu+6Cyy8/79ubNIHp03XhTdOmOuee6oxUyk1bvb+DIlUD83cJcZuKNANe5PsVG0gYrc9F2v8b4cxG7Kn44Ycf/g71pUuXMnz4cI8rCoyLL76YQ4cOsXTpUi677DKvy/GfXLlOh/iqVfDJJ7rF0nXXwc8/n/etTZtquG/cqJ/v2XPGF5PvQHxrQAsmrj6Vpe/2zJ4rLzkrNwNnERKS0rMxqr8/Qnkz6z59+gggcXFxcuLECa/LCZjKlSsLIN9//73XpQTevn0ijz8ukiePbrrdqZPIgQPnfcusWSJxcSLVqons2ZP85PHD8tML1wggzzbNLvJ4fpFJPUWOHw783yEEdR61QDqPWuB1GVkK6dzM2n7cJjt06BDOOUaNGsWzzz7L0aNHiYmJ8bqsgLjyyitZtWoVc+bMoUGDBl6XE3h588Jjj+nSl0ce0ZudUvZn3bcv1bc0bw7TpsG6ddCsmd4bdSwpissHTwfggUa57W5PE7rSk/7+/gi1EfuyZcsEEEA2bNjgdTkB1aZNGwFk+vTp6X5PxI3MkpL0zz17RPLmFbnpJpH161N96cyZItmzi1x+uQjkF0BOvH6lyC/fiIxqKLJychALN1kdNmJPvxw5ctChw3UkJiaG/wXE8/j++++ZMWMGEydOpF27dul6z2XtbmNinwYs/nVv5Cz1TFndFBUFvXrpHHylStCzpy6LOUPLltpg8uefTwJzWLhwPTH9voPSDaH3N97eiWvMOdidp1lIYmIiO3bsoHjx4mm+9sCBA+TNmxeA3DXbc0GzXtQtfQEQgasf/vwTnn8e3ngDTp3SBe1nfI82b95MuXJ3Eh39BdWqRTN7NlxwgYf1mizLWgqYTPvoo4+46aabAFi3bh2Pzt8LhGagp/wG4ZfaduzQO5Vuv10fjxsHTZqQVKwY3333HUeONKJDB90HZPZsyJ/f91MakxFBbSngnGvlnFvvnNvknHvAH8c0wXfq1CmKFy/OTTfdRKNGjUhKSqJChQpel3VOCaMXsvjXvf6bJipS5HSo//UX9O0L5coRNWAAjcqVo3VrnbVZuVJvbj3HdVdjPOdzsDvnooHXgdbApcCNzrlLfT2uCa6FCxcSExPDtm3bmDNnDvPnz8clz0VP6F0/5EbrIsKP7z3D78+34+T/ArDXXYECsHYt9OgBo0dD2bJw9920qbOHKVN0OXzz5hbuJjT5Y8ReB9gkIr+IyAlgPHCtH45rgkBE6NChAw0aNCBnzpwcP36cpk2bel3Wec2dO5eoqCh+XfA5FVvezJU1Kgfmh0/JkhrqGzZA167w1ltw/Djt2sHUKcLPP9vI3YQmfwR7MWDLGY+3Jj/3D865Xs65Jc65Jbuty1JI+O2334iKimLatGmMHj2aQ4cO/d2LPRTt37+f2NhYmjVrRt68eTl48CBVr+sb+BOXLq1tgrdsgWL6T7vdmPZMueZdli8XWraE/fsDX4Yx6eWPYE+tM9a/rsiKyBgRqSUitQoVKuSH0xpf7Nu3j9KlSwPae71Xr14eV3R+gwYNIl++fJw8eZJvv/2Wffv2kStXruBOE6VcLT12DPLk4ZpPbmVydBd+WnKKlk1PWbibkOGPYN8KXHzG4+JAACY9jT/lzZuX77//HhEJ6d7rq1evxjnH8OHD6devHyLCVVdd5W1RcXHw4YewahXtr3VMSurI0qVCqysOcOCAt6UZA/4J9h+BS5xzpZ1zsUAX4DM/HNcEkHMuLNoJ7N+/n6JFi7J3715ef/11r8v5p0svhfHjuXbFU0xq/AZL1uemVSs4sHA1lvDGSz4Hu4icAu4CvgLWAhNFJLR3ejVho0GDBmzbto38obxovEoVOnw9gAkTHD/+KLRucpyDparAs8/CoUNeV5cpb775JvvsqnDY8ss6dhGZISLlRaSsiDztj2MaE26uv+JYPE4AABoWSURBVB7Gj3csPnk5rZnBwYee0Quvw4bB4cNel3duZ2wc8sorr+Cco1evXqxYscLrykwmRVSvmIjpZWLCVseOMGGCY9GBy2hTbTsHq18Fgwdrk/dQtXMF4x5KwEVFcffdd1OnTh2OHDlCw4YNva7MZFLE7KCUMHohi37ZAyL+vc3cmAzq2BHGj4cuXXLTtsFUZsxbQq6rkjf+GD1ad3Xq1Usvwnpp11qYcjtvzt1Er2nHKFsghqUPXU7ebu/qdoMmbEXMiH3b8m/444X2HFg8xetSjKFTJ/joI1iwANo+VotDR6P1C7Nn6/6s5crByJEa8l7JVxIKlqdTueP8NTg3m/rnIm+Jy/R547mkpCR+OavbaHpFRLCPHTuW70c9SO7CpWiWcHtI3gJvzi1Sp9A6d9ZVkf/9L7RtmzzNPnkyzJ2rc+933gmXXKKNx7wQGw+NBpM/TrggV3bbOCSE/Pzzz0RHR1O2bNlMvT8igr1jx46sWLGC1o9/RFR0xMwuZQkJoxfy9Zef8dW4VyIy4BMSUgn3Jk3g229h1iy9kzVPHn3xgQNw8mRwC9y1BopUh25T9M9dtqDNE2dcwO7ZsyfVq1cHdGe3TB7PdlAy3ti0adPfO1flqt468nZqOsNHH4lERYlcfbXIoUNnfCFlNycRkX79RMqWFXn3XZGTJ4Neo/HQ9uWy7Z5cf///MGzYsFRfRjp3ULLhrQm648ePU7t2bVauXAlAq8c/Jk/hkhE9fXbjjToou/lmuOYa+PxziI/n9G5OAG3a6KR8jx7wzDO6T2tCAkRHe1W2CbTkC9gc3k2xl3R0vv3ZWhS5pa1Ph42IqRgTPkSEuLg4Vq5cyXvvvYeIkKdw1rhY17UrvPcefPONhvuRI2e9oG1bWLYMpk7VFTM33QQPPeRJrSZIki9gc3gPmwfkQh7PT5FyVXy+gG07KJmgW7x4MXXq1Pm733tW88EHcMstOtX+2WfJI/ezJSXBlClQo4b2gl+5UtsHX3ed7tVqIseutTCyHkTHQuIJ6LcYLqyY6kuDuoOSMRlRt27dLBvqAN266a57X38N114LR4+m8qKoKLjhBg110P1YO3XSoJ82Ted1TGQIwAVsG7Eb45Fx46BnT2jaVEfu570nKDFR73p64gnYuFED/plnoGXLoNVrvGcjdmNCXPfuun/H3LnnGbmniI7WOfc1a+Ddd+F//9M1lKCjdxvBmzNYsBvjoR494O23Yc4cnT4/diyNN2TLpj8R1q+HBx/U56ZPhyuv1INYwBss2I3xXM+eup3qV19Bhw7pCHeAmJjTV10TE+GPP3R37UaNYP78QJZrwoAFuzEh4NZbT4f79denM9xTXHcdbNoEI0bA5s3QuLEuuzFZlgW7MSHittvgzTfhyy+1Q2SG+oNlz669ZzZvhuHDoVUrff74cVi8OCD1mtBlwW5MCLn9dhgzBmbM0JF7hps/xsVp98iuXfXxO+9AvXp6V+uPP/q9XhOaLNhDQCg2vwrFmrKKO+7Qtu0zZmRi5H62bt3g+efhhx+gTh295XXZMr/VakKTT+vYnXPDgGuAE8BmoKeIpLlRoq1jP63TiPl89+235ChTk7qlLwC83yCk0+vfMvv9V8ld8xquqFo+JGrKikaPhj59NIsnT4bYWB8OdvAgvPoqvPii9oJfvPiffWpMWAjWOvbZQGURqQpsAB708XhZyltvvcWU/o3ZNfkJr0v525AhQ5hyV6PkDUvsf3wv9e6te3FMn643oZ444cPBcueGoUPh11/h/fc11Pfs0RH9qlV+q9mEBp+6O4rIrDMeLgI6+VZO1lK0aDHuvPNOdlXuinMuJEbFhw8f5q677uLPy24MmZqysr59dWn6nXfqxh0TJ/o4cs+bVz9Ap2Q++0y3ekpI0G6SFVPvUWLCiz/n2G8FvjzXF51zvZxzS5xzS3bv3u3H04avNm1aM2LEiJDqmzJixAhee+21kKopq+vXT1cyTpum+evTyP1MLVroCP6BB/TXgssu077Cp0756QTGK2nOsTvn5gCFU/nSUBGZlvyaoUAt4HpJx6S9zbEbkz579uxh3rx53HDDDbz2GgwYoMvWJ0zQe5T8ZvduGDYMduzQqRo9ORQs6MeTGF+ld449zakYEWmWxom6A+2ApukJdWNM+pw4cYJChQoB2se+f3+dlhk4ELp00Z5gfgv3QoXghRdOtyRYvx6qVdMR/NChUKqUn05kgsGnqRjnXCtgCNBeRM7eNsAYk0kiQvbs2YF/7ns5YAC8/LLuxXHjjQHYIjVlCi5fPr16+957UL68Ls/ZssXPJzOB4usc+wggNzDbObfcOTfKDzUZk+UVLVoUgN9//52cOXP+42t33w0vvaT7cHTtGqD9ry+6CF55Re9kvf12bUNZuTJkdnNlE1S+roop569CjDGqY8eO7Ny5kwULFlCiRIlUXzNokM6a3HuvDrI/+kgbP/pd8eK65nLIEF37niuXPv/qq7oGs0iRAJzU+MruPDUmRCSMXkiVa3szdepU3n//ferXP/9S03vugf/8ByZN0lbtAV3MUrKkrrcEnX+/5x4oU0Z/suzaFcATm8ywYDcmBCSMXsjiX/ey6rMxVGhxE9MPl03X++69VxezTJyo9xoFZaVihQqwbp0G/fDhULq0jugPHAjCyU16WLAbE0JK3D+NatffmaH33HefLmiZMCGIy9DLldO9/das0Sby48ad3mTbFsd5LhCzcsaYDJrQu/7fTdcyc7fv/fdrng4Zovk6blyA5tzPVqECfPihjtZz5dKfKvXrazfJQYN0dY0JOhuxGxMhBg+GZ5/VC6ndu+vGSkGTJ4/+eeCAzsc/+aSuff+//7MpGg/41N0xs+zOU2MC59ln4aGH9ILquHG6D3bQLV8Ojz+ufRDy59eNty+91INCIovf7jw1xoSXBx/UaZmhQ3Va5p13PAj36tXh009h6VLdrbtCBX1+8WJdD3/W2nzjXxbsxkSghx6CpCR45BFd5z52rEcj95o19QPg6FGde8+WTS8G9O0LOXJ4UFTkszl2YyLUww/rVPd77+l+qkGdc09NjhzaRbJKFV2nWaaM3uiUoZ27TXpYsJuQU//2J7nhjf96XUZEeOQRneoeN047A3ge7g0awJw58M032vt94ECdrjF+ZVMxJqQkjF7I4nHP8Ed0URJG69yBbfbhm8ce0zn3J57QaZm33jq95NwzDRvCvHm62UeNGvrck09qj5qePX3cTcR4/Z/XmH+5+O7xxOS3HiT+9Pjj8OijeiH1jjt0/j0kpIR6UhLMnatdJMuX1wuuAeluljVYsJuQMqF3feqVu4i6pS9gQu/6Nlr3o8cf13n3sWO1I2/IhDvorxDz58MXX2hv+Ntv16ma777zurKwZMFuTBbhnM52DB2q0zF9+oRYuDunq2Z++EH3Yi1QAJLbF/PXXyFwgSB82By7CTk2Sg8c5/RmUBF45hkdKI8cGQJz7mdyDq65Rj9S9OgBmzbprx033BBiBYce++4Yk8U4B089pTcyjR4Nd94ZYiP31PTooQvxu3TRLfumTAmDor1jwR5CEkYv/LsRlDGB5Bw8/bTeJzRqFNx1V4g3ZezYEVas0I1eT52CTp3gtde8ripk2VRMiOg8agGz3n+NXNVbkzBan7MpCRNIzmlfGRFt+wvw+uuntz0NOVFRkJCgoT5+vM7Hg66JP3RIH4ds8cHllxG7c+4+55w45wr643hZkXOO2MLliM5pbU59JgJrPgvxIWhocA6ee07b/r7xBvTvHwbftuho7XCWP78+fvllaNcO6tWDmTPD4C8QeD6P2J1zFwPNgT98LyfrmtC7PglnfG58sHMFTLwZen8HRap6XU3Icw6ef17z8D//0cevvhpGg99Jk7Rvwv/9H7Rurf3gn38errrK68o8448R+8vAYMB+TBpPDXr1Q3rUv5AT4zpR7tVDuKLVcM7R9Kp6rFmzxuvyQppzOh1z770wYoTe6R82A9+YGG2Gs2GDXjDYskW37oOse4FVRDL9AbQHXkn+/Deg4Hle2wtYAiwpUaKEGONPnUctkFzl6wkgJx/NK//tGS83Vo6R7DHRAsgHH3zgdYlhISlJ5J57REBk4EB9HHaOHRM5cUI/f+01kauvFvnuO29r8hNgiaQjm9OcinHOzQEKp/KlocBDQIt0/gAZA4wB3WgjPe8xJiNyXNaES2rUJZt7hzol4riiRDbotxgurOh1aWHDOZ2OSUrSfaqdg5deCqNpGYDs2U9/nisXrF2r0zLNm2vDnPqRP9WZ5lSMiDQTkcpnfwC/AKWBn51zvwHFgWXOudR+CBgTUBN616dxy3bcUrcwFKlOTPdPoEh12LXa69LCTkqYDxyo4X7vvWE0LXO2Hj3gl1/0p9Xy5dpd8v77va4q4DJ98VREVgIXpjxODvdaIrLHD3UZkykLczTi7t4P6IPe33hbTBhzThebiOifKSP5sBq5p4iP159OffroBYRayTvL7dkDf/xxuhFZBLF17CZi2Goi/3JOR+wiOoKPitILrGEZ7qDb8Q0Zcvrx8OF6l1aHDtqqoFo1z0rzN78Fu4iU8texjDGhwTl45ZV/LoV8/vkwDvcz3X+/9n1/6SXdn7VjRw34ypW9rsxnNmI3xpxXyrr2pCQYNkxH7s8+GwHhnjevNqnv31/D/ZVX9PnJk72tyw+sV4wxJk3O6fR03746Yn/ooTC+oHq2/Pn15qZff9ULCqDr4G++WdfGhyELdmNMuqSEe58+2oZg6NAICnfQ/u8XX6yfL18OU6dCpUq6smbzZk9LyygLdmNMukVFaaOwXr10OubhhyMs3FN06aLLJO++GyZMgAoVtAVmmLBgN8ZkSFSUNgy74w7drOPRRyM03C+6CF58UQP+rrv0ZifQv+zOnd7Wlga7eGqMybCoKG3LIqKbdjinN3WG/QXV1BQpoksjU8ybp83G7rhDdyspVsy72s7BRuzGmEyJitIdmG67Ta89PvGE1xUFSfny0LOn/uXLltXpmhAbwVuwG2MyLSoKxoyBW2/VYM8S4V68uP66smGD9oUfMUJ7wYfQZtsW7MYYn0RFwZtv6iD28cfhySe9rihISpeGt9/WpZGjRukGIKdO6VXlPd52VrFgN8b4LCoK3npLVwY+9phOzWQZ5cpBq1b6+YIFug60dGn9c+9eT0qyYDfG+EVKuN9yi66UefppryvyQMOGsHo1tG2rS4ZKldKfdMeOBbUMC3ZjjN9ER8PYsXrT5sMPa7ZlOZUq6WbbK1ZAixZ6o1NMjH4tSDs62XJHY4xfRUfDO+/oUsihQwGSqFlzNi1btvS6tOCqUkX7zhw5ot+UAwegZk290nzXXZA7d8BObSN2Y0LIyZMn2bJli9dl+Cw6Gt59F7p0SWTo0ChatZpHUlbdfzQ+Xv/cvx8qVtRGO6VLaw/kw4cDckoLdmNCxKlTp4iNjaVkyZJel+IXhw7tZ/z4WOBD4DmGDcvicXPxxTB9OixeDLVra2/4MmVg1y6/nyqLf6eNCQ2JiYnEJM/D/vnnnx5X47vt27eTL18+IIm9e9tz443wwAPa9jfLq1MHvvwSvv8ebr8dLkzeiG7ePL9dZLVgN8Zjnd/4nmzZ9HLXzp07KVSokMcV+Wb9+vUUS77N/vjx4+TPn5v33tO+WoMH64YdBt1/NWXp0M6d0LKl3sn6+utw/LhPh7ZgN8ZDCaMXMn/uLACueW4aAz79xeOKfLNo0SIqVqwI6G8hsbGxAGTLBu+/DwkJunHRiy96WWUIKlwYvvpKp2buugsuuURv6T1xIlOH8znYnXP9nXPrnXOrnXMv+Ho8Y7KaHGVqUmLwdHLkC++R+ueff079+vW56KKLEBGiov4ZL9mywQcfQOfOcN99p/e0MMkaN4Zvv4VZs7SxWL9+sHVrpg7l03JH51xj4Fqgqogcd85d6MvxjMlqJvSuT8LohX9/Hq7Gjh3LbbfdRoMGDfj+++/P+bps2eDDD3Up5D336HODBgWpyHDgHDRvDs2a6Y1OZcro83fdpXPz6eTrOva+wHMichxARPx/edcYE7ISRi9kzYx3WPXZm9x000188MEHab4nJdyTkjTcndMGieYMzp3eVPvwYVi4EBYtSvfbfQ328sBVzrmngWPAfSLyo4/HNCZLCdeResLohcwa/xb75o+lYsubOXlV33S/NyYGPv5YL6gOGqTtCAYMCGCx4SxnTliyBHbv1s0/0iHNYHfOzQEKp/Klocnvzw/UA2oDE51zZUT+vZ+Kc64X0AugRIkS6SrOGBPaonPl54Lmfal63c0Zfm9MjN55n5AAAwdquIfR7nPB5dzpZZHpeXkqGZyBc7mZ6FTM/OTHm4F6IrL7fO+rVauWLFmyJNPnNcaEBn9cHzhxQi+oTpumrc3vvNNf1UUe59xSEamV1ut8XRXzKdAk+YTlgVjA20bExpiwEhsLEyfCtdfqiH3kSK8rCn++zrGPBcY651YBJ4DuqU3DGGMik7+uD6SE+w036IjdOeib/il7cxafgl1ETgDd/FSLMSYLi42FSZOgUyddwu0c9OnjdVXhye48NcaEjJRwb9dOR+yjR3tdUXiyYDcR6/pX5nLdy7O8LsNkUPbs2sa8bVsdsY8Z43VF4cc22jARKWH0Qma+9iCJR/aRED8OCN/14llR9uwwZQpcfz307q3TMnfc4XVV4cOC3USs/I1vJTrXBV6XYTLpzHDv1UvD/fbbva4qPFiwm4g0oXd9Es743ISnuDjdMvS663TE7hzcdpvXVYU+m2M3xoS0uDj45BNtV37HHbqfqjk/G7GbiGUj9cgRFweffqo3Md12m47ce/TwuqrQZSN2Y0xYSAn3Zs3g1lth3DivKwpdFuzGmLCRI4f2lGnaFHr21F2ZzL9ZsBtjwkpKuDdpAt27665M5p8s2I0xYSc+Hj77THeT695dN+4wp1mwG2PCUnw8TJ8OjRrBLbfARx95XVHosGA3xoStlHBv2BBuvll3ZTIW7MaYMJczJ3z+OVx1FXTrprsyZXUW7MaYsJczJ3zxBVx5Jdx0E0yY4HVF3rJgN8ZEhJRwv+IKDfeJE72uyDsW7MaYiJErF8yYAQ0aQNeu2ts9K7JgN8ZElJRwr18fbrxRe7tnNRbsxpiIkxLu9epBly7a/jcr8SnYnXPVnXOLnHPLnXNLnHN1/FWYMcb4Indu+PJLqFtXw33qVK8rCh5fR+wvAE+ISHXg0eTHxhgTElLCvXZtSEjQ9r9Zga/BLkCe5M/zAtt9PJ4xxvhVnjwwc6aGe+fO2iEy0jkRyfybnasEfAU49IdEAxH5/Ryv7QX0Sn5YGViV6RMHT0Fgj9dFpIPV6T/hUCNYnf4WLnVWEJHcab0ozWB3zs0BCqfypaFAU+AbEZninOsM9BKRZmme1LklIlIrrdd5zer0r3CoMxxqBKvT3yKtzjR3UDpfUDvn3gMGJj+cBLyV7gqNMcYEhK9z7NuBRsmfNwE2+ng8Y4wxPvJ1z9M7gFecc9mAY5yeQ0/LGB/PGyxWp3+FQ53hUCNYnf4WUXX6dPHUGGNM6LE7T40xJsJYsBtjTITxLNjDqR2Bc66/c269c261cy5k7651zt3nnBPnXEGva0mNc26Yc26dc26Fc+4T51w+r2s6k3OuVfJ/503OuQe8ric1zrmLnXPznHNrk/89Dkz7Xd5wzkU7535yzn3udS3n4pzL55ybnPzvcq1zrr7XNaXGOTco+b/3Kufcx865uPO93ssRe1i0I3DONQauBaqKyGXAfzwuKVXOuYuB5sAfXtdyHrOByiJSFdgAPOhxPX9zzkUDrwOtgUuBG51zl3pbVapOAfeKSCWgHnBniNYJuhR6rddFpOEVYKaIVASqEYL1OueKAQOAWiJSGYgGupzvPV4Ge7i0I+gLPCcixwFEZJfH9ZzLy8Bg9PsakkRkloicSn64CCjuZT1nqQNsEpFfROQEMB79gR5SRGSHiCxL/vwgGkTFvK3q35xzxYG2hPC9Lc65PEBD4G0AETkhIvu8reqcsgE5klcgxpNGXnoZ7HcDw5xzW9BRcMiM3s5SHrjKObfYOfeNc6621wWdzTnXHtgmIj97XUsG3Ap86XURZygGbDnj8VZCMDDP5JwrBVwOLPa2klQNRwcaSV4Xch5lgN3AO8lTRm8553J6XdTZRGQbmpF/ADuA/SIy63zv8XUd+3mlox3BoDPaEbwNpNmOIBDSqDMbkB/9tbc2MNE5V0aCvE40jRofAloEs55zOV+dIjIt+TVD0SmFD4NZWxpcKs+F7G8/zrlcwBTgbhE54HU9Z3LOtQN2ichS59zVXtdzHtmAGkB/EVnsnHsFeAB4xNuy/sk5lx/97bE0sA+Y5JzrJiIfnOs9AQ32cGlHkEadfYGpyUH+g3MuCW0YtDtY9cG5a3TOVUH/g//snAOd3ljmnKsjIjuDWCJw/u8lgHOuO9AOaBrsH45p2ApcfMbj4oTo9KBzLgYN9Q9FJBS7jF8BtHfOtQHigDzOuQ9EpJvHdZ1tK7BVRFJ+45mMBnuoaQb8KiK7AZxzU4EGwDmD3cupmHBpR/ApWh/OufJALCHUBU5EVorIhSJSSkRKof9Ya3gR6mlxzrUChgDtReSI1/Wc5UfgEudcaedcLHpx6jOPa/oXpz+93wbWishLXteTGhF5UESKJ/977AJ8HYKhTvL/I1uccxWSn2oKrPGwpHP5A6jnnItP/u/flDQu8gZ0xJ6GzLYjCLaxwFjn3CrgBNA9xEaa4WQEkB2YnfzbxSIR6eNtSUpETjnn7kLbUEcDY0VktcdlpeYK4GZgpXNuefJzD4nIDA9rCmf9gQ+Tf5j/AvT0uJ5/SZ4mmgwsQ6cwfyKN1gLWUsAYYyKM3XlqjDERxoLdGGMijAW7McZEGAt2Y4yJMBbsxhgTYSzYjTEmwliwG2NMhPl/cJPfqLMAV/0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "lr = LogisticRegression(solver='lbfgs')\n",
    "lr.fit(Xtest.numpy(), ytest.numpy())\n",
    "beta_train = a_tch.detach().numpy().flatten()\n",
    "beta_test = lr.coef_.flatten()\n",
    "b_train = b_tch.squeeze().detach().numpy()\n",
    "b_test = lr.intercept_[0]\n",
    "hyperplane = lambda x, beta, b: - (b + beta[0] * x) / beta[1]\n",
    "\n",
    "Xtrain_np = Xtrain.detach().numpy()\n",
    "Xtrain_grad_np = Xtrain_grad.numpy()\n",
    "ytrain_np = ytrain.numpy().astype(np.bool)\n",
    "\n",
    "plt.figure()\n",
    "plt.scatter(Xtrain_np[ytrain_np, 0], Xtrain_np[ytrain_np, 1], s=25, marker='+')\n",
    "plt.scatter(Xtrain_np[~ytrain_np, 0], Xtrain_np[~ytrain_np, 1], s=25, marker='*')\n",
    "\n",
    "for i in range(m):\n",
    "    plt.arrow(Xtrain_np[i, 0], Xtrain_np[i, 1],\n",
    "              Xtrain_grad_np[i, 0], Xtrain_grad_np[i, 1])\n",
    "\n",
    "plt.xlim(-8, 8)\n",
    "plt.ylim(-8, 8)\n",
    "\n",
    "plt.plot(np.linspace(-8, 8, 100),\n",
    "         [hyperplane(x, beta_train, b_train)\n",
    "          for x in np.linspace(-8, 8, 100)], '--', color='red', label='train')\n",
    "plt.plot(np.linspace(-8, 8, 100),\n",
    "         [hyperplane(x, beta_test, b_test)\n",
    "         for x in np.linspace(-8, 8, 100)], '-', color='blue', label='test')\n",
    "plt.legend()\n",
    "plt.savefig(\"data_poisoning.pdf\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "@webio": {
   "lastCommId": null,
   "lastKernelId": null
  },
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
